* =============================================================================
* 波士顿房价多元回归预测分析
* 文件名: boston.do
* 数据源: dataset/boston.csv
* 作者: 自动生成
* 日期: 2025/11/16
* =============================================================================

clear all
set more off
set seed 12345

* =============================================================================
* 1. 数据导入与预处理
* =============================================================================

* 导入CSV数据
import delimited using "dataset/boston.csv", clear

* 显示数据结构
describe

* 显示前10行数据
list in 1/10

* 检查缺失值
misstable summarize

* 变量标签
label variable crim "犯罪率"
label variable zn "住宅用地比例(>25000平方英尺)"
label variable indus "工业用地比例"
label variable chas "是否临河(1=是,0=否)"
label variable nox "氮氧化物浓度(ppm)"
label variable rm "平均房间数"
label variable age "1940年前建造房屋比例"
label variable dis "到就业中心加权距离"
label variable rad "高速公路可达性指数"
label variable tax "财产税率(每万美元)"
label variable ptratio "师生比"
label variable b "黑人比例"
label variable lstat "低收入人口比例"
label variable medv "房价中位数(千美元)"

* =============================================================================
* 2. 描述性统计分析
* =============================================================================

* 基本统计描述
summarize

* 详细统计描述
summarize, detail

* 分类变量频数分析
tabulate chas

* =============================================================================
* 3. 数据可视化
* =============================================================================

* 设置图形输出格式
graph set window fontface "Arial"
graph set eps fontface "Arial"

* 房价分布直方图
histogram medv, frequency title("房价中位数分布") ///
    xtitle("房价中位数(千美元)") ytitle("频数")

* 房价箱线图
graph box medv, title("房价中位数箱线图") ///
    ytitle("房价中位数(千美元)")

* 散点图矩阵（主要变量与房价）
graph matrix medv rm lstat dis nox, half title("主要变量与房价关系")

* =============================================================================
* 4. 相关性分析
* =============================================================================

* 计算相关系数矩阵
pwcorr medv rm lstat dis nox crim age tax ptratio, star(0.05)

* 相关性热力图（需要安装heatmap包）
* 如果没有安装，可以使用以下命令安装：
* ssc install heatmap, replace
* heatmap medv rm lstat dis nox crim age tax ptratio, title("变量相关性热力图")

* =============================================================================
* 5. 多元线性回归模型
* =============================================================================

* 基础多元回归模型（包含所有变量）
regress medv crim zn indus chas nox rm age dis rad tax ptratio b lstat

* 逐步回归（向前选择）
stepwise, pr(0.1): regress medv crim zn indus chas nox rm age dis rad tax ptratio b lstat

* 逐步回归（向后剔除）
stepwise, pe(0.1): regress medv crim zn indus chas nox rm age dis rad tax ptratio b lstat

* 基于逐步回归结果建立最终模型
* 假设逐步回归选择了以下重要变量
regress medv rm lstat dis nox crim ptratio rad tax b

* =============================================================================
* 6. 模型诊断
* =============================================================================

* 残差分析
predict residuals, residuals
predict fitted_values, xb
predict std_resid, rstandard

* 残差正态性检验
swilk residuals

* 残差vs拟合值图
scatter residuals fitted_values, yline(0) title("残差vs拟合值") ///
    xtitle("拟合值") ytitle("残差")

* Q-Q图
qnorm residuals, title("残差Q-Q图")

* 杠杆值和Cook距离
predict leverage, leverage
predict cooksd, cooksd

* 识别异常值（|标准化残差| > 2）
gen outlier = abs(std_resid) > 2
list if outlier == 1

* =============================================================================
* 7. 多重共线性检验
* =============================================================================

* 计算方差膨胀因子（VIF）
* 需要安装estatokg包
* ssc install estatokg, replace
estat vif

* =============================================================================
* 8. 模型改进与比较
* =============================================================================

* 对数变换模型
gen log_medv = log(medv)
regress log_medv rm lstat dis nox crim ptratio rad tax b

* 二次项模型
gen rm_sq = rm^2
gen lstat_sq = lstat^2
regress medv rm lstat rm_sq lstat_sq dis nox crim

* 交互项模型
gen rm_lstat = rm * lstat
regress medv rm lstat rm_lstat dis nox crim

* 模型比较（AIC、BIC）
estimates store model1
regress log_medv rm lstat dis nox crim ptratio rad tax b
estimates store model2
estimates stats model1 model2

* =============================================================================
* 9. 预测与验证
* =============================================================================

* 数据分割：80%训练集，20%测试集
set seed 12345
gen random = runiform()
sort random
gen train = _n <= _N * 0.8

* 在训练集上建立模型
regress medv rm lstat dis nox crim ptratio rad tax b if train == 1
estimates store train_model

* 在测试集上进行预测
predict pred_test if train == 0, xb

* 计算预测性能指标
gen actual = medv if train == 0
gen error = actual - pred_test
gen abs_error = abs(error)
gen sq_error = error^2

* 计算MAE、RMSE、MAPE
summarize abs_error if train == 0
local mae = r(mean)
display "平均绝对误差(MAE): " `mae'

summarize sq_error if train == 0
local rmse = sqrt(r(mean))
display "均方根误差(RMSE): " `rmse'

gen ape = abs_error / actual * 100 if train == 0
summarize ape if train == 0
local mape = r(mean)
display "平均绝对百分比误差(MAPE): " `mape' "%"

* =============================================================================
* 10. 结果可视化
* =============================================================================

* 实际值vs预测值散点图
scatter actual pred_test if train == 0, title("实际值vs预测值") ///
    xtitle("预测值") ytitle("实际值")

* 残差分析图
scatter pred_test error if train == 0, yline(0) title("预测值vs残差") ///
    xtitle("预测值") ytitle("残差")

* 变量重要性图（标准化系数）
regress medv rm lstat dis nox crim ptratio rad tax b
* 计算标准化系数
foreach var in rm lstat dis nox crim ptratio rad tax b {
    summarize `var'
    local sd_`var' = r(sd)
}
summarize medv
local sd_medv = r(sd)

* 显示标准化系数
display "=== 标准化系数 ==="
foreach var in rm lstat dis nox crim ptratio rad tax b {
    local beta = _b[`var']
    local std_beta = `beta' * `sd_`var'' / `sd_medv'
    display "`var': " `std_beta'
}

* 手动创建标准化系数图
matrix beta = e(b)'
matrix list beta

* =============================================================================
* 11. 保存结果
* =============================================================================

* 保存预测结果
preserve
keep if train == 0
keep actual pred_test error abs_error
order actual pred_test error abs_error
export delimited using "output/boston_predictions.csv", replace
restore

* 保存模型结果
estimates save "output/boston_model.ster", replace

* 生成分析报告
log using "output/boston_analysis.log", replace
display "=== 波士顿房价多元回归分析报告 ==="
display "数据集大小: " _N " 个观测值"
display "模型R²: " e(r2)
display "调整R²: " e(r2_a)
display "F统计量: " e(F)
display "模型显著性: " e(p)
display "平均绝对误差(MAE): " `mae'
display "均方根误差(RMSE): " `rmse'
display "平均绝对百分比误差(MAPE): " `mape' "%"
log close

* =============================================================================
* 12. 清理临时变量
* =============================================================================

drop random train outlier residuals fitted_values std_resid leverage cooksd
drop log_medv rm_sq lstat_sq rm_lstat pred_test actual error abs_error sq_error ape

display "=== 分析完成 ==="
display "结果文件已保存到 output/ 目录"
display "包括预测结果、模型文件和分析日志"

* =============================================================================
* 程序结束
* =============================================================================
